/*!
 * @file exponential.h
 * @ingroup MathUtil
 *
 * Routines for generating random numbers.
 *
 */



#ifndef _QT_EXPONENTIAL_H_
#define _QT_EXPONENTIAL_H_

#include <cstdlib>
#include <cmath>
#include <cassert>

#include <boost/random.hpp>

#include "mathutil.h"
#include "randomutil.h"
#include "uniform.h"

namespace qt {

namespace random {

/*!
 *  This class provides functions and default generators
 *  associated with the exponential distribution.
*/
class Exponential {

public:

    /*!
     * Cumulative distribution.
     *
     * @param x is a real number
     * @returns Prob(X<x)..
     *
     * History:
     *
     */
    static double cdf(double x,double lambda) {
        assert(x>0);
        return 1-std::exp(-1*lambda*x );
    }


    /*!
     * Inverse of the cumulative distribution function.
     *
     * @param u a real number strictly between 0.0 and 1.0
     * @returns x the largest number such that P(X<x) < u.
     *
     * History:
     *
     */
    static double invCdf(double u,double lambda) {
        return -1.0/lambda * std::log(1-u);
    }


    /*!
     * The density function of the  distribution with respect to dx.
     *
     * @param x a real number
     * @returns Prob( x-dx < X < x+dx )/dx
     *
     * History:
     *
     */
    static double pdf(double x, double lambda) {
        assert(x>0);
        return lambda*std::exp(-1*lambda*x);
    }

public:

    Exponential(double lambda = 1.0) : _lambda(lambda) {
        assert(lambda>0);
    }


    /*!
    * Cumulative distribution.
    *
    * @param x is a real number
    * @returns Prob(X<x)..
    *
    * History:
    *
    */
    double cdf(double x) {
        return cdf(x,_lambda);
    }


    /*!
     * Inverse of the cumulative distribution function.
     *
     * @param u a real number strictly between 0.0 and 1.0
     * @returns x the largest number such that P(X<x) < u.
     *
     * History:
     *
     */
    double invCdf(double u) {
        return invCdf(u,_lambda);
    }


    /*!
     * The density function of the  distribution with respect to dx.
     *
     * @param x a real number
     * @returns Prob( x-dx < X < x+dx )/dx
     *
     * History:
     *
     */
    double pdf(double x) {
        return pdf(x,_lambda);
    }



    template <typename GeneratorType>
    class GeneratorWrapper {

    public:

        typedef DistributionGenerator<GeneratorType,
                Exponential> Generator;

        typedef WrappedGenerator
        < GeneratorType,
        DistributionGenerator<GeneratorType, Exponential > >
        BoundGenerator;


    };


    typedef
    GeneratorWrapper<qt::random::DefaultGeneratorType>::Generator
    DefaultGenerator;

    typedef
    GeneratorWrapper<qt::random::DefaultGeneratorType>::BoundGenerator
    DefaultBoundGenerator;


    static
    DefaultGenerator getDefaultGenerator(double lambda) {
        DefaultGeneratorType & dgen =
            DefaultGeneratorType::getDefaultGenerator();
        Exponential e(lambda);
        DefaultGenerator gen(dgen,e);
        return gen;
    }

    template <typename GeneratorType>
    static
    typename GeneratorWrapper<GeneratorType>::Generator
    getGenerator(GeneratorType & ugen, double lambda) {
        Exponential e(lambda);
        typename GeneratorWrapper<GeneratorType>::Generator gen(ugen,e);
        return gen;
    }

    /*! Returns an instance of a generator that owns
     *  its own base generator.  Reseeding this generator
     *  will not affect other generators.
    */
    template<typename GeneratorType>
    static
    typename GeneratorWrapper<GeneratorType>::BoundGenerator
    getBoundGenerator(double lambda,
                      unsigned long seed = default_seed) {

        Exponential e(lambda);

        typename GeneratorWrapper<GeneratorType>::BoundGenerator wgen(e);

        return wgen;
    }


    double get_lambda() {
        return _lambda;
    }

    void set_lambda(double lambda) {
        _lambda = lambda;
    }

private:
    double _lambda;

};

}; // namespace qt::random

}; // namespace qt::


#endif
